from scipy.special import gdtrix, ndtri
# import matplotlib.pylab as plt
import matplotlib.pyplot as plt
import numpy as np
import sqlite3
import pandas as pd
from sqlalchemy import create_engine

# 解决matplotlib显示中文乱码问题
plt.rcParams['font.family'] = ['sans-serif']
plt.rcParams['font.sans-serif'] = ['SimHei']  # windows, Linux
# plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']  # Mac

# 数据库文件位置


# 单站日降水量数据读取，安溪县
df = pd.read_excel("/home/matrix/workspace/PycharmProjects/Disasters/rain_and_flood/SURF_ANXI_DAY_1961_2019.xls")
df = df.sort_values(by='PRE_24h', ascending=False, axis=0).reset_index(drop=True)
mean_pre = df.iloc[:, 1].mean()
df.fillna(value=mean_pre, inplace=True)
Q = df.iloc[:, 1].astype("int")


def xtrans(plist):
    """
    x坐标轴转换成海森机率格纸格式
    :param plist: plist传入的是一个p值列表数据，plist必须是概率值，一定是大于1的数字，
    如6%，取6，返回的xplist也是一个列表数据
    :return xplist: 符合海森机率格式的x坐标轴
    """
    xzero = ndtri(0.0001)
    xplist = []
    for i in range(len(plist)):
        xplisti = ndtri(plist[i] / 100)
        xplisti -= xzero
        xplist.append(xplisti)
    return xplist


def ljxs(plist, Cs):
    """
    求离均系数fp
    :param plist: plist传入的是一个p值列表数据，plist必须是概率值，一定是大于1的数字，
    如6%，取6，返回的xplist也是一个列表数据
    :param Cs:
    :return:
    """
    aa = 4 / Cs ** 2  # b值，反映的是形状（shape），在水文计算中为α（alpha），不要与scipy中的官方文档中的a值混淆
    tp = []  # 为求离均系数fp，tp是过渡
    fp = []
    for i in range(len(plist)):  # 离均系数fp
        tpi = round(gdtrix(1, aa, 1 - plist[i] / 100), 3)  # 官方文档中的a值为scale，本次采用了标准伽马分布，所以取a=1
        fpi = round(Cs * tpi / 2 - 2 / Cs, 3)
        tp.append(tpi)
        fp.append(fpi)
    return fp  # 返回的fp也是一个列表数据


# 坐标轴的绘制
pstandard = [0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0, 20.0, 25.0,
             30.0, 40.0, 50.0, 60.0, 70.0, 75.0, 80.0, 90.0, 95.0, 97.0, 99.0, 99.5, 99.9]
x_axis = [str(i) + '%' for i in pstandard]
x = xtrans(pstandard)
y_axis = np.linspace(0, round(max(Q) * 1.2, 2), 10, dtype=np.float32).round(1)
print(x_axis, y_axis, sep='\n')
print(x)
# 建立海森图纸张，需要知道Q的最大值来确定横坐标的上限
fig = plt.figure(figsize=(12, 6))
ax1 = fig.add_subplot(111)
for i in range(len(x)):
    plt.vlines(x[i], 0, round(max(Q) * 1.2, 2), 'g', '--')
for j in range(len(y_axis)):
    plt.hlines(y_axis[j], 0, max(y_axis), 'g', '--')
plt.xticks(x, x_axis, color='blue', rotation=90)
plt.yticks(y_axis, y_axis, color='blue', rotation=0)
plt.ylim(0, round(max(Q) * 1.2, 2))
plt.xlim(0, round(max(x) + 0.1, 2))
plt.tick_params(labelsize=10)
ax1.set_title("安溪县1961-2019年日降水量PIII分布适线图")
ax1.set_xlabel("概率设定值P%")
ax1.set_ylabel("降水量（0.1mm/d）")

# 绘制原始数据的散点图
qj = np.mean(Q)
qcv = np.std(Q) * np.sqrt(len(Q) / (len(Q) - 1)) / qj  # 无偏估计的Cv值,标准差除均值
# Q.sort(reverse=True)
print("降水量Q（0.1mm）:\n {}".format(Q))
pq = pd.Series([(i + 1) * 100 / (len(Q) + 1) for i in range(len(Q))])  # Q的经验频率 P=m/(n+1)
print("Q的经验频率 P=m/(n+1):\n{}".format(pq.astype("int")))
pqx = xtrans(pq)
plt.scatter(pqx, Q)
# 绘制预测的概率曲线，需要知道Cs,Cv,qj（均值），一般通过倍比n=Cs/Cv来计算Cs=Cv*n
n = 2  # 是自定义值
Cs = n * qcv
qpredict = [qj * (qcv * m + 1) for m in ljxs(pstandard, Cs)]
plt.plot(x, qpredict, 'r-', lw=3)
plt.text(6.5, max(Q) * 1.1, r'mean_value={:.2f},Cs={:.2f}, Cv={:.2f}'.format(qj, Cs, qcv), size=15, \
         family="Monospace", color="black", weight="normal", verticalalignment='top', \
         horizontalalignment='right', bbox=dict(facecolor="w", alpha=1))
plt.show()
